/*
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#include <mpi.h>

// -- matriz.h
double** alocar_matriz(int altura, int largura);
void desalocar_matriz(double **M, int altura);
double** ler_matriz(FILE *file, int intervalSize, int intervalStart, int *nread, int *J_ORDER);
double* ler_b(FILE *file, int *nread, int *J_ORDER, int intervalSize, int intervalStart);
void imprimir_matriz(double **M, int altura, int largura);
// -- main.c
int main(int argc, char **argv)
{
    int J_ORDER = 0, J_ROW_TEST = 0, J_ITE_MAX = 0; //Ordem da matriz, linha a ser testada, maximo de iteracoes
    double **A, *B, *D, *newD, *aux; //Matriz A, vetor B, vetor resultado D, e auxiliar para D
    double J_ERROR, ctrl, error_ctrl = 1, max = 0, soma = 0, test = 0, max_diff = 0;  //variaveis de controle para o algoritmo de Jacobi
    int i, j, k;         //variaveis de iteracao
    double dtest =  0;    //variavel para testar o resultado em J_ROW_TEST
    int nread = 0;       //numero de linhas lidas
    int myrank, nprocs;  //rank deste processo, numero de processos
    int intervalSize;    //tamanho do intervalo que cada processo ira calcular (linhas da matriz)
    int intervalStart;   //quantas linhas pular ateh começar a primeira linha do intervalo de cada processo
    int intervalEnd;     //fim do intervalo.
    int *sendcounts, *recvcounts;     //ler explicacao na documantacao do MPI, funcao MPI_alltoallv(...);
    int *sdispls, *rdispls;        //ler explicacao na documantacao do MPI, funcao MPI_alltoallv(...);
    int numthreads;
    FILE *file;

    MPI_Init(&argc, &argv);                 //inicializa MPI
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank); //pega o rank deste processo
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs); //numero total de processos

    file = fopen(argv[1], "r");
    if(file == NULL){
        printf("Impossivel abrir o arquivo.\n");
        exit(-1);
    }

    if(argc != 3){
        printf("usage: jacobiP <filename> <num_threads>");
        exit(-1);
    }

    sscanf(argv[2], "%d", &numthreads);
    if(numthreads < 1){
        printf("Numero de threads invalido\n");
    	exit(-1);
    }

    if(nprocs == 0){
        printf("Erro na criacao dos processos\n");
    }

    fscanf(file,"%d\n",&J_ORDER);            //leitura de variaveis
    fscanf(file,"%d\n",&J_ROW_TEST);
    fscanf(file,"%lf\n",&J_ERROR);
    fscanf(file,"%d\n",&J_ITE_MAX);

    if(nprocs > J_ORDER){
        if(myrank == 0)
            printf("\nEste programa não suporta o uso de um numero de processos maior que o a ordem da matriz.\n");
        exit(-1);
    }

    intervalSize = J_ORDER/nprocs;           //determinacao dos intervalos da matriz que cada processo ira calcular
    intervalStart = intervalSize * myrank;
    intervalEnd = intervalStart + intervalSize;

    recvcounts = (int*) malloc (nprocs * sizeof(int));  //alocacao de vetor de quantidades
    rdispls = (int*) malloc (nprocs * sizeof(int));     //alocacao de vetor de offsets
    sendcounts = (int*) malloc (nprocs * sizeof(int));  //alocacao de vetor de quantidades
    sdispls = (int*) malloc (nprocs * sizeof(int));     //alocacao de vetor de offsets

    if(recvcounts == NULL || rdispls == NULL || sdispls == NULL || sendcounts == NULL){
        printf("Falha ao alocar memoria\n");
        exit(-1);
    }

    for(i = 0; i < nprocs ; i++){
        recvcounts[i] = intervalSize;    //seta quantidades a serem lidas de cada no quando MPI_alltoallv() for chamado
        rdispls[i] = i * intervalSize;   //seta offset a partir do qual escrever quando MPI_alltoallv() for chamado
        sendcounts[i] = intervalSize;    //seta quantidades a serem enviadas a cada no
        sdispls[i] = 0;                  //seta offset do endereco inicial (de D) a ser enviado ao no de rank i
    }

    recvcounts[nprocs-1] += J_ORDER % nprocs;  //ajusta a quantidade do ultimo nó para ler mais o resto da divisao, se houver

    if(myrank == nprocs-1) {              //se este processo for o ultimo, processa o trecho dele + o resto
        intervalSize += J_ORDER % nprocs;
        for(i = 0; i < nprocs ; i++)      //como o ultimo processo calcula o trecho dele + o resto, ele precisa enviar
            sendcounts[i] = intervalSize; //uma quantidade maior de elementos para os demais a cada iteracao
    }

    A = ler_matriz(file, intervalSize, intervalStart, &nread, &J_ORDER); //Matriz principal
    B = ler_b(file, &nread, &J_ORDER, intervalSize, intervalStart);      //Vetor B
    D = (double*) calloc(J_ORDER , sizeof(double));     //D vetor com os valores encontrados na iteracao atual (inicialmente 0,0,0...0)
    newD = (double*) calloc(J_ORDER , sizeof(double));     //auxiliar para calcular o erro.
    aux = (double*) calloc(intervalSize , sizeof(double));   //auxiliar para o D da proxima iteracao

    if(D == NULL || aux == NULL || newD == NULL){
        printf("Falha ao alocar memoria\n");
        exit(-1);
    }

    omp_set_num_threads(numthreads);

#pragma omp parallel shared(A,B,ctrl) private(k,i)
    {
    
    #pragma omp for schedule(static)
    for (k = 0; k < intervalSize; k++) // "divide a matriz pela diagonal" e zera a diagonal
    {
    	ctrl = A[k][k + intervalStart]; //diagonal eh deslocada para a direita intervalStart posicoes (0 no no mestre)
    	if(ctrl == 0){
            printf("Matriz invalida\n");
            exit(-1);
        }
    	B[k] /= ctrl;                    //divisao dos componentes de B pelo respectivo elemento da diagonal (B torna-se B*)
    	for (i = 0; i <J_ORDER; i++)
			A[k][i] = A[k][i]/ctrl;   //divisao dos componentes da linha de A pelo respectivo elemento da diagonal (A torna-se A*)

        A[k][k + intervalStart] = 0;   //zera a diagonal, mais eficiente do que uma comparacao dentro do loop
    }

    }   //fim da regiao paralela (pragma)

    for (k = 0; k < J_ITE_MAX && error_ctrl; k++) //enquanto k < numero de iteracoes maximo e erro nao satisfeito
    {	
	#pragma omp parallel shared(A,D,aux) private(j,i,soma,test)
	{          
            #pragma omp for schedule(static)   
	for (i = 0; i < intervalSize; i++){     //algoritmo de jacobi
	    soma = test = 0;
	    for (j = 0; j < J_ORDER; j++){
            soma += A[i][j] * D[j];
            test += fabs(A[i][j]);
            }

            aux[i] = B[i] - soma;  //uso de vetor auxiliar para nao alterar o D enquanto ainda esta sendo usado
            if (test >= 1){        //aux agora contem um trecho dos valores de D
               printf("Sistema divergente\n");
               exit(-1);
            }
	
	}         
	#pragma omp barrier
	}//fim da regiao paralela	




        /* Funcao principal do programa, faz a troca de todos as partes de D entre todos os nós:
        a partir de aux, cada no envia sendcounts[J] variaveis do tipo double
        e recebe recvcounts[J] variaveis do tipo double do processo de rank J na posicao
        newD[rdispls[J]] no contexto MPI_COMM_WORLD */

        MPI_Alltoallv(aux, sendcounts, sdispls, MPI_DOUBLE, newD, recvcounts, rdispls, MPI_DOUBLE, MPI_COMM_WORLD);

        /*note que esta funcao envia de todos para todos, inclusive o proprio processo, assim newD ja sera atualizado
        com o valor de aux pois cada processo enviara aux para todos inclusive para si proprio, e salvara o que
        receber em newD*/

        /*newD agora contem os valores atualizados por todos os nos do vetor D*/
        max = 0;
        max_diff = 0;
        for (i = 0; i < J_ORDER; i++){

            if(fabs(newD[i]) > max)
                max = fabs(newD[i]);  //max pega o maior entre os valores de newD[i] (valores parciais do vetor resposta B atualmente)
            if(fabs(newD[i]) - fabs(D[i]) > max_diff)
                max_diff = fabs(newD[i]) - fabs(D[i]);  //pega a maior das diferencas entre cada newD[i] e D[i]
            D[i] = newD[i];                             //atualiza D com o valor de newD para a proxima iteracao
        }

        if(max != 0 && max_diff / max < J_ERROR) //checa se satisfaz o erro
            error_ctrl = 0;                      //se sim zera error_ctrl que quebra o loop

        if(myrank == 0)
            if(k % 100 == 0)
                printf("Iteracao %d alcancada (maximo: %d)\n", k, J_ITE_MAX);

	}

    MPI_Barrier(MPI_COMM_WORLD);             //garante que todos os processos chegaram ateh aqui

    /*apenas o no que processou a linha pedida ira imprimir o resultado*/

    if(J_ROW_TEST >= intervalStart && J_ROW_TEST < intervalEnd){
        for (i = 0; i < intervalSize; i++) // volta a diagonal para 1 (estava com 0 anteriormente para o algoritmo de jacobi)
            A[i][i + intervalStart] = 1;

        for(i = 0; i < J_ORDER; i++){      // dtest recebe o valor da equacao da linha J_ROW_TEST avaliada com as variaveis encontradas
            int temp = (J_ROW_TEST % intervalSize);
            dtest += D[i] * A[temp][i];
        }

        printf("\n---------------------------------------------------------\n");
        printf("Iterations: %d\n", k);
        printf("RowTest:%d => %lf =? %lf\n", J_ROW_TEST, dtest, B[J_ROW_TEST % intervalSize]);
        printf("---------------------------------------------------------\n");
    }

    //Limpando
    desalocar_matriz(A, intervalSize);
    if(B != NULL)
        free(B);
    if(newD != NULL)
        free(newD);
    if(D != NULL)
        free(D);
    if(aux != NULL)
        free(aux);
    if(sendcounts != NULL)
        free(sendcounts);
    if(recvcounts != NULL)
        free(recvcounts);
    if(sdispls != NULL)
        free(sdispls);
    if(rdispls != NULL)
        free(rdispls);

    fclose(file);
    MPI_Finalize();
    return 0;
}

//Aloca matrizes
double** alocar_matriz(int altura, int largura)
{
    int i;
    double** M = (double**) malloc(altura * sizeof(double*));
    for(i = 0; i < altura; i++){
        M[i] = (double*) malloc(largura * sizeof(double));
        if(M[i] == NULL){
            printf("Erro na alocacao da matriz\n");
            exit(-1);
        }
    }
    if (M == NULL)
        exit(-1);
    return M;
}

//Desaloca matrizes
void desalocar_matriz(double **M, int altura)
{
    int i;
    for(i = 0; i < altura; i++)
        if(M[i] != NULL)
            free(M[i]);
    if (M != NULL)
        free(M);
}

void imprimir_matriz(double **M, int altura, int largura)
{
    int i, j;
    for(i = 0; i < altura; i++){
        for(j = 0; j < largura; j++)
            printf("%lf ", M[i][j]);
        printf("\n");
    }
}

//Leitura e alocação de matrizes (chama alocar_matriz)
double** ler_matriz(FILE *file, int intervalSize, int intervalStart, int *nread, int *J_ORDER)
{
    double **M;
    int i, j;
    char c;

    while(*nread < intervalStart){    //a partir da 4a linha, le \n ateh chegar no comeco do intervalo
        while((c = fgetc(file)) != '\n');
        (*nread)++;
    }

    M = alocar_matriz(intervalSize, *J_ORDER);  //aloca uma matriz de intervalSize x ordem da matriz

    for (i = 0; i < intervalSize; i++){       //faz a leitura do trecho da matriz
        for(j = 0; j < *J_ORDER; j++)
            fscanf(file,"%lf",&M[i][j]);

    }
    *nread = intervalStart + intervalSize;   //posiciona nread no final do trecho lido, para coincidir com o file pointer
    return M;
}


//Leitura do vetor B
double* ler_b(FILE *file, int *nread, int *J_ORDER, int intervalSize, int intervalStart)
{
    int i;
    double *B = (double*) malloc (intervalSize * sizeof(double));
    char c;

    while(*nread <= *J_ORDER + intervalStart){    //a partir de onde ler_matriz deixou o file pointer, le \n ateh chegar no comeco do intervalo
        while((c = fgetc(file)) != '\n');
        (*nread)++;
    }

    for (i = 0; i < intervalSize; i++){     //faz a leitura do vetor B (apenas o trecho necessario)
        fscanf(file,"%lf\n",&B[i]);
    }

    return B;
}
